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•  We  model  lithiation  in  crystalline  Si  anode  based  on  discrete  phase  boundary  motion. 

•  Anisotropic  phase  boundary  motion  is  taken  into  account. 

•  Microscopic  model  is  derived  to  explain  size-dependent  fracture  during  lithiation. 

•  We  estimate  the  work  of  fracture  of  lithiated  Si. 
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Silicon  (Si)  nanostructures  are  attractive  candidates  for  electrodes  for  Li-ion  batteries  because  they 
provide  both  large  specific  charging  capacity  and  less  constraint  on  the  volume  changes  that  occur  during 
Li  charging.  Recent  experiments  show  that  crystalline  Si  anodes  expand  highly  anisotropically  through 
the  motion  of  a  sharp  phase  boundary  between  the  crystalline  Si  core  and  the  lithiated  amorphous  Si 
shell.  Here,  we  present  a  microscopic  model  to  describe  the  size-dependent  fracture  of  crystalline  Si 
nanopillars  (NPs)  during  lithiation.  We  derive  a  traction-separation  law  based  on  the  plastic  growth  of 
voids,  which,  in  turn,  is  used  in  a  cohesive  zone-finite  element  model.  The  model  allows  for  both  the 
initiation  of  cracking  and  crack  growth.  The  initial  size  and  spacing  of  the  nanovoids,  assumed  to  be 
responsible  for  the  fracture,  together  with  the  computed  facture  toughness,  are  chosen  to  conform  to 
recent  experiments  which  showed  the  critical  diameter  of  Si  NPs  to  be  300-400  nm.  The  anisotropy  of 
the  expansion  is  taken  into  account  and  that  leads  naturally  to  the  observed  anisotropy  of  fracture.  The 
computed  work  of  fracture  shows  good  agreement  with  recent  experimental  results  and  it  may  be 
possible  to  use  it  to  describe  the  failure  for  other  loading  and  geometries. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

High-capacity  lithium-ion  batteries  (LIBs)  have  attracted  much 
attention  as  a  key  element  for  portable  electronic  devices  and 
electric  vehicles  1,2].  The  search  for  anode  materials  with  high 
charging  capacity  has  identified  silicon  (Si)  as  one  of  the  most 
promising  candidates,  because  it  has  an  exceptionally  high  specific 
charging  capacity  of  4200  mAh  g-1  [3].  However,  the  associated 
huge  volume  expansion  during  lithiation,  can  cause  pulverization 
and  capacity  loss  [4,5],  thus  impeding  the  development  of  Si  an¬ 
odes.  To  avoid  the  mechanical  degradation  due  to  lithiation- 
induced  stresses,  various  nanostructures  have  been  suggested  as 
possible  electrodes  for  improved  cycling  performance  and  fracture 
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resistance.  Examples  include  nanowires  [6],  thin  films  [7,8],  core¬ 
shell  structures  [9],  and  hollow  structures  [10,11]. 

Recent  experiments  have  shown  that  lithiation  of  crystalline  Si 
differs  from  that  for  amorphous  Si:  the  lithiation  process  for  crys¬ 
talline  Si  involves  a  two-phase  reaction,  with  an  atomistically  sharp 
phase  boundary  existing  between  the  amorphous  LixSi  shell  and 
the  pristine  crystalline  Si  core  [12,13].  Moreover,  the  phase 
boundary  moves  faster  in  the  <110>  direction  of  crystalline  Si  than 
other  directions  [12,14,15].  These  phenomena  differ  from  the  lith¬ 
iation  of  amorphous  Si,  where  lithiation  is  governed  mainly  by 
diffusion. 

In  an  effort  to  understand  the  lithiation  process  for  crystalline  Si, 
and  especially  to  understand  the  stress  evolution  and  fracture, 
several  authors  have  analyzed  the  stress  evolution  and  deformation 
associated  with  interface  reaction-controlled  lithiation.  Zhao  et  al. 
[16]  developed  an  analytical  model  for  concurrent  interface 
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reaction  and  plasticity  in  nano-spheres.  Pharr  et  al.  [17]  developed  a 
kinetics-based  model  for  anisotropic  motion  of  phase  boundaries 
by  specifying  phase  boundary  motion  and  computing  morphology 
changes  associated  with  lithiation,  similar  to  those  found  experi¬ 
mentally  [18].  By  adjusting  the  diffusivity  of  lithium  ions  to  simu¬ 
late  a  sharp  boundary,  Yang  et  al.  [19]  also  developed  a  model  to 
show  the  anisotropic  shape  changes.  In  addition,  McDowell  et  al. 
[20]  have  suggested  a  kinetics  model  to  account  for  the  effect  of 
stress  on  the  interface  reaction  rate. 

In  the  present  study,  we  have  developed  a  microscopic  model  to 
explore  size-dependent  fracture  of  crystalline  Si  nanopillars  (NPs) 
during  lithiation.  As  a  first  approximation,  we  model  the  phase 
boundary  motion  by  specifying  a  temperature  field,  using  the 
analogy  between  thermal  expansion  and  lithiation-induced 
swelling.  Using  this  model,  we  can  simulate  the  anisotropic 
expansion  and  accompanying  stress  evolution  during  the  lithiation 
of  crystalline  Si  NPs  with  various  crystallographic  orientations. 

To  model  the  crack  initiation  and  propagation,  we  first  derive  a 
traction-separation  law  based  on  a  microscopic  description  of  the 
plastic  growth  of  voids.  From  the  traction-separation  law,  the  work 
of  fracture  (Gc)  is  estimated  by  taking  the  initial  size  and  spacing  of 
the  nanovoids  such  that  a  critical  diameter  of  Si  NPs  of 300-400  nm 
is  predicted,  consistent  with  experiment. 

We  also  performed  experiments  to  determine  the  critical 
diameter  for  fracture  of  Si  NPs  during  lithiation  for  different  crys¬ 
tallographic  orientations.  We  used  the  scanning  electron  micro¬ 
scope  (SEM)  to  observe  the  shapes  of  individual  pillars  after 
lithiation,  as  done  previously  18].  These  experiments  are  then 
directly  compared  with  the  modeling  results  and  the  comparison 
leads  to  a  good  agreement  on  the  size  effect  on  fracture.  This 
characterization  method  determines  a  critical  diameter  for  fracture 
and  also  serves  as  a  standard  to  judge  the  validity  of  our  modeling. 

2.  Stress  evolution  during  lithiation 

2.1.  Isotropic  expansion  model 

Recent  transmission  electron  microscopy  (TEM)  studies  showed 
that  a  sharp  phase  boundary  exists  between  crystalline  Si  and 
amorphous  LixSi,  thus  creating  a  core-shell  structure  during  lith¬ 
iation  of  a  crystalline  Si  NP  [13,21,22].  As  a  first  approximation,  we 
start  with  the  isotropic  expansion  model.  To  derive  an  approximate 
analytical  solution  for  the  stresses  during  lithiation,  we  divide  the 
domain  between  the  crystalline  Si  core  and  the  amorphous  LixSi 
shell,  as  shown  in  Fig.  1A.  In  this  model,  crystalline  Si  and  the 
amorphous  LixSi  are  assumed  to  be  elastic  and  elastic-plastic, 
respectively.  When  the  phase  boundary  moves  during  lithiation, 


(C) 


Fig.  1.  Schematic  of  lithiation  of  a  crystalline  Si  nanopillar  (A)  At  a  given  time,  pristine 
silicon  core  and  LixSi  alloy  shell.  The  thin  layer  (blue  color)  stands  for  the  atomistically 
sharp  phase  boundary  only  where  lithiation  occur  at  the  given  time  duration.  (B)  Stress 
state  of  the  lithiated  Si  alloy  shell.  (C)  Stress  state  of  pristine  Si  core.  (For  interpretation 
of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version 
of  this  article.) 


the  interface  reaction  involving  the  breaking  of  Si-Si  bonds  and  the 
forming  of  Li-Si  bonds  is  assumed  to  occur  only  in  an  atomistically 
thin  layer  (blue  color  in  Fig.  1  A).  Due  to  the  volume  expansion  of  the 
thin  layer  between  the  core  and  shell,  the  surrounding  LixSi  shell 
can  be  considered  to  be  in  the  same  stress  state  as  a  pressurized 
cylindrical  tube,  with  the  core  being  subjected  to  a  uniform  radial 
pressure,  as  shown  in  Fig.  IB  and  C.  [16,19,23].  Following  Hill’s 
treatment  of  a  pressurized  tube  with  the  Tresca  yield  criterion  [24], 
the  stresses  in  the  LixSi  shell  would  then  be 

(t)m  =  -'"(;)•  (tL,’1  -'"(?)•  <lab> 

after  assuming  that  the  tube  deforms  fully  plastically  due  to  the 
large  volume  expansion.  on ,  oqq  and  Y  are  radial,  tangential  and 
yield  stress,  respectively,  b  is  the  outer  radius  of  tube  and  r  is  the 
distance  from  the  center.  The  corresponding  stresses  in  the  Si  core 
would  be 


where  a  is  the  radius  of  the  Si  core.  From  this  analytical  solution  for 
the  hoop  stress,  we  can  see  that  a  tensile  stress  should  develop  at 
the  surface. 

To  account  for  the  stress  state  in  the  interfacial  layer,  we  use  a 
finite  element  package,  ABAQUS  (2010  version),  for  lithiation- 
induced  swelling  with  plasticity.  In  this  model,  we  adopt  an 
isotropic  elastic  and  perfectly  plastic  model,  using  the  mechanical 
properties  shown  in  Table  1.  Using  the  analogy  between  diffusion- 
induced  swelling  and  thermal  expansion,  we  model  the  phase 
boundary  by  specifying  the  temperature  field.  The  total  strain  rate 
at  any  point  is  composed  of  the  elastic  (£e),  and  plastic  (£p)  strain 
rates  as  well  as  the  transformation  strain  rate  (e1)  due  to  lithiation. 
Elastic  properties  are  assumed  to  be  isotropic,  obeying  Hooke’s  law, 
and  the  conventional  J-2  flow  rule,  without  hardening,  is  used  for 
plastic  strain  calculations.  The  amorphous  LixSi  shell  expands  iso¬ 
tropically  according  to  the  specified  temperature  distribution  with 
a  linear  transformation  strain  of  0.4,  which  leads  to  about  a  300% 
volume  expansion.  The  reaction  front  thickness  is  taken  to  be  much 
smaller  than  the  radius.  For  the  out-of  plane  deformation,  we  as¬ 
sume  a  plane  strain  condition. 

With  the  numerical  model  described  above,  stresses  have  been 
computed  for  the  lithiation  of  a  Si  NP  with  a  diameter  of  100  nm. 
The  hoop  stresses  are  plotted  in  Fig.  2.  It  shows  tensile  stress 
developed  easily  at  the  surface  due  to  huge  volume  expansion  at 
the  interfacial  boundary,  which  provides  a  clear  explanation  for 
why  cracking  is  initiated  at  the  surface  of  the  NPs  during  lithiation 
of  crystalline  Si  [25],  rather  than  from  the  center. 

2.2.  Anisotropic  expansion  model 

Our  recent  experimental  study  revealed  that  lithiation  of  crys¬ 
talline  Si  occurs  faster  in  the  <110>  direction  than  other  directions 
[18,26].  To  model  the  anisotropic  phase  boundary  motion,  we 


Table  1 

Material  properties  and  operating  parameters. 


Material 

Description 

Symbol  [dimension] 

Value 

References 

Crystalline  Si 

Young’s  modulus 

E  [GPa] 

185 

[36] 

Yield  point 

y  [GPa] 

7 

[37] 

Poisson’s  ratio 

V 

0.28 

[36] 

LixSi  alloy 

Young’s  modulus 

E  [GPa] 

35 

[38] 

Yield  point 

y  [GPa] 

1 

[39] 

Poisson’s  ratio 

V 

0.22 

[40] 
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Fig.  2.  Hoop  stress  from  the  isotropic  expansion  model  of  a  crystalline  Si  pillar  with 
initial  radius  of  50  nm.  This  plot  shows  the  stress  across  the  cross  section  from  the 
center  of  the  pillar  to  the  surface  at  five  different  times  during  the  lithiation  process; 
legend  shows  the  degree  of  lithiation. 


adjust  the  temperature  field  according  to  different  crystallographic 
orientations  of  the  Si  core.  Following  experimental  observations, 
the  phase  boundary  on  the  <110>  plane  is  assumed  to  move  5 
times  faster  than  that  on  the  <100>  plane  and  the  velocities  in  all 
other  directions  are  assumed  to  be  much  smaller. 

With  the  prescribed  temperature  fields,  we  compute  stresses 
and  morphology  changes  in  Si  NPs  with  <100>,<110>,  and  <111  > 
axial  orientations.  We  observe  that  morphology  changes  according 
to  the  prescribed  temperature  fields  are  in  a  good  agreement  with 
experimental  observations  [18],  as  shown  in  Fig.  3.  As  expected,  the 
stress  concentration  occurs  midway  between  adjacent  <110> 
planes,  which  coincides  with  the  fracture  sites  in  experiments. 
Detailed  morphologies  and  stress  evolution  are  presented 
in  Supplementary  movies  1,  2,  and  3. 

Supplementary  video  related  to  this  article  can  be  found  at 
http://dx.doi.org/10.1016/jjpowsour.2013.12.137. 


3.  Microscopic  model  for  fracture 

Based  on  the  TEM  observation  that  nanometer-sized  pores  can 
form  during  amorphization  of  Si  [27  ,  we  developed  a  microscopic 
description  of  fracture  based  on  the  plastic  growth  of  nanovoids.  To 
explore  size-dependent  fracture  of  crystalline  Si  NPs  upon  lith¬ 
iation,  we  use  cohesive  elements  in  the  finite  element  method. 
Dugdale  [28]  and  Barenblatt  [29]  first  introduced  the  concept  of 
cohesive  zones  near  the  crack  tip  and  Needleman  30]  developed  a 
finite  element  method  for  the  implementation  of  that  concept.  In 
the  present  model,  spherical  nanovoids  are  assumed  to  form  in 
amorphous  LixSi  at  a  certain  stress.  During  lithiation,  the  nanovoids 
grow  and  eventually  link  up  and  lead  to  fracture,  as  shown  in 
Fig.  4A,  B.  For  simplicity,  the  interaction  between  the  voids  is 
ignored  and  the  spherical  voids  are  assumed  to  grow  in  response  to 
the  hydrostatic  tensile  stress.  As  shown  in  Fig.  4C,  we  model  the 
cohesive  zone  by  considering  hypothetical  bars  subjected  to 


♦  t  t 


(A)  (B)  (C) 

Fig.  4.  Schematic  of  the  plastic  void  growth  near  the  crack  tip.  (A)  Nanovoids  are 
formed  with  the  radius  of  a  and  the  spacing  of  2/?  under  the  loading  and  (B)  eventually 
linked  up  and  lead  to  failure.  (C)  Hypothetical  bar  to  model  the  hydrostatic  stress  state 
near  the  void.  The  total  volume  of  the  bar  is  conserved,  because  the  volume  increase  at 
the  edges  is  same  as  the  volume  decrease  due  to  the  void  growth.  a0  is  the  initial  radius 
of  void  and  w,  wD  are  the  separation  distance  at  current  and  at  the  void  formation, 
respectively. 


Fig.  3.  Hoop  stress  evolution  from  the  anisotropic  expansion  model  of  a  crystalline  Si  pillar  with  different  crystallographic  orientation  (A)  <100>,  (B)  <110>,  and  (C)  <111  >.  The 
initial  radius  of  pillars  is  50  nm.  Blue  triangles  indicate  the  fracture  sites  in  the  experiments.  Values  below  figures  denote  the  degree  of  lithiation  and  all  figures  are  plotted  with  the 
same  contour  scale.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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uniaxial  tensile  deformation,  where  the  void  is  under  a  hydrostatic 
stress  state,  so  that  it  would  grow  isotropically,  keeping  the 
spherical  shape.  Afterward,  when  the  radius  of  the  plastic  domain  is 
same  as  the  void  spacing,  the  stiffness  of  the  cohesive  elements 
would  start  to  be  degraded. 


3.1.  Plastic  void  growth 


Following  the  analysis  of  McClintock  [31]  for  a  cylindrical 
void,  we  assume  that  spherical  voids  start  to  form  when  the 
hydrostatic  stress  reaches  a  critical  value,  which  is  later 
computed  from  the  assumed  void  radius  a  and  void  spacing  2/?. 
We  start  by  deriving  the  stress  field  of  a  growing  spherical  void 
in  an  ideally  plastic  solid.  Due  to  the  spherical  symmetry,  all 
variables  depend  only  on  the  radial  direction.  The  equilibrium 
equation  and  kinematic  relations  in  the  spherical  coordinate 
system  can  be  expressed  by 


dan  2  (ffrr  -  a99)  _  dii  .  dii 

—j - 1 - —  U,  £rr  —  ~r~ ,  £0$  —  — , 

dr  r  dr  00  r 


(3a,b,c) 


where  u  is  the  displacement  in  the  radial  direction  in  the  spherical 
coordinate  system. 

From  J- 2  plasticity,  strain  rates  are  expressed  by 
f'i  =  \  (fff  -  ffm)  Y  (i  =  rr’ 99 -  #)  (4) 


where  om  is  the  mean  stress  and  i  is  the  effective  strain  rate. 

Using  the  strain  relations  and  equilibrium  condition,  the  radial 
stress  can  be  expressed  by 


&rr(P) 


e 

4Y  fe99-errdr 

3  J  i  r’ 

a 


(6) 


after  using  the  traction  boundary  condition  at  the  inner  surface. 
From  the  kinematic  relation  and  incompressibility  condition,  the 
strain  rate  can  be  expressed  by 

err  =  -2a  ^  (7a) 

'£ee  =  =  ^^3  (7b) 


(7c) 


where  a  is  the  rate  of  change  of  the  void  radius.  Substituting  these 
expressions  into  Equation  (6),  we  get 

M£)  =  2Ylog(£)+2j,  (8) 

where  ys  is  the  surface  energy.  The  last  term  is  added  to  account 
for  the  effect  of  the  surface  energy  of  the  void  surface.  This  rela¬ 
tion  gives  the  hydrostatic  tension  stress  needed  to  maintain 
plastic  growth  of  the  void;  below  we  assume  that  this  also  cor¬ 
responds  to  the  tensile  stress  in  the  hypothetical  tensile  bar,  as  if 
an  incompressible  fluid  were  present  above  and  below  the 
growing  void, 


3.2.  Traction- separation  law 

Prior  to  the  void  formation,  the  stress-strain  relations  for  hy¬ 
pothetical  bars  subjected  to  uniaxial  tensile  deformation,  are  as 
described  in  Appendix  A.  To  obtain  the  traction-separation  law  for 
this  regime,  we  adopt  these  relations  and  estimate  the  total  strain 
in  z-direction  as  the  opening  displacement  divided  by  the  void 
spacing  (/?),  as  follows  : 

a  =  (2/i  +  A)etotal  =  (2/i  +  A)^  (Elastic  regime)  (9a) 

a  =  ^  +  A^total+|y 

=  &-  +  a)  j  +  (Plastic  regime)  (9b) 

where  a  is  the  traction  and  w  is  the  separation  displacement  of  the 
cohesive  zone  and  A,  g  are  Lame’s  constants.  By  equating  these 
relations,  we  can  find  the  traction  (op)  and  the  separation  distance 
(wP)  corresponding  to  the  onset  of  the  plasticity,  as  follows  : 

op  =  ^1  Wp  =  (Onset  of  plastic  deformation) 

(10a,b) 

We  find  the  traction  (gd)  at  the  point  of  void  formation  by 
setting  a  =  a0,  the  initial  void  size,  in  Equation  (8),  and  find  the 
separation  distance  (wD)  at  this  point  by  equating  the  tractions  in 
Equations  (8)  and  (9b): 


aD  =  2Y\og(— ^  +2— ,  wD  =  — 3-/2  (Void  formation) 

V*o  J  3A  +  2m 

(11a, b) 

From  the  volume  conservation  during  the  void  growth,  and 
again  imagining  an  incompressible  fluid  above  and  below  the  void, 
as  shown  Fig.  4C,  we  can  find 

(w- wd)-tt/32  =  |rc(a3  -  a$)-  D2) 

Substituting  Equation  12  into  Equation  (8),  the  traction-sepa¬ 
ration  law  during  plastic  void  growth  is  then 


3(w  -  wD) 

4^ 


3(w  -  wD)  (jf\ 
4/3  VaV 


(Void  growth  regime). 


(13) 


As  we  have  noted,  the  hydrostatic  tension  stress  can  be  equated 
to  the  traction  of  the  cohesive  zone  if  the  hypothetical  bar  is 
assumed  to  behave  like  a  fluid,  meaning  that  uniaxial  loading  is 
transmitted  as  the  hydrostatic  stress  around  the  void.  Since  a 
cohesive  element  completely  loses  its  load-carrying  capacity  when 
the  void  size  is  same  as  the  spacing  of  voids,  there  is  a  maximum 
separation  distance  (Wf)  which  can  be  computed  from  Equation 
(12)  by  setting  a  =  P,  as  follows: 


4  / 

wf  =  -  y  +wD  (Maximum  separation)  (14) 

3  V  P  J 

Finally,  the  traction-separation  law  we  will  use  is  plotted  in 
Fig.  5  for  various  values  of  the  initial  void  size  and  their  spacing.  In 
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Fig.  5.  Traction-separation  law  computed  from  the  microscopic  model  of  plastic  void  growth.  (A)  Traction-separation  law  and  (B)  accompanying  fracture  toughness  with  constant 
spacing  of  3  nm  and  variable  void  radius.  (C)  Traction-separation  law  and  (D)  accompanying  fracture  toughness  with  constant  initial  void  radius  of  1.5  nm  and  variable  void  spacing. 


the  traction-separation  law,  at  a  given  void  spacing  of  3  nm,  Gc 
increases  with  decreasing  void  size  (in  Fig.  5A,  B).  In  addition,  Gc 
increases  with  increasing  void  spacing  (in  Fig.  5C,  D)  at  a  fixed 
initial  void  radius  of  1.5  nm.  These  could  be  understood  by  the  fact 
that  more  supporting  material  exists  in  the  cohesive  zone  with 
decreasing  void  size  and  increasing  void  spacing. 

3.3.  Size-dependent  fracture  modeling  with  cohesive  zone  elements 

With  the  prescribed  traction-separation  law,  we  explore  the 
size-dependence  of  fracture  in  <100>  oriented  Si  NPs  during 
lithiation.  Cohesive  elements  with  zero  thickness  are  implanted 
along  45°,  the  tilted  direction  between  two  (110)  planes,  where  a 
crack  is  expected  to  propagate  during  lithiation  due  to  the  stress 
concentration  in  this  plane.  Satisfying  the  equilibrium  conditions, 
cohesive  elements  convey  traction  to  the  adjacent  elements,  based 
on  the  traction-separation  law,  and  can  simulate  crack  initiation  as 
well  as  its  propagation.  A  detailed  explanation  of  how  to  simulate 
crack  propagation  with  cohesive  elements  is  given  in  Appendix  B. 

With  cohesive  elements,  we  performed  simulations  to  estimate 
the  Gc  which  would  lead  to  the  critical  size  of  300-400  nm,  which 
was  experimentally  observed  [18]  in  <100>  oriented  Si  NPs  during 
lithiation.  By  varying  the  initial  void  size  and  their  spacing,  we 
could  determine  whether  cohesive  elements  completely  fail  at  the 
surface  according  to  the  prescribed  traction-separation  law.  For 
instance,  with  very  low  Gc,  we  could  see  the  cohesive  element  at 
the  surface  fails  and  the  crack  propagates  to  the  center,  but  for  this 
Gc,  fracture  occurs  easily  even  for  small  Si  NPs  less  than  100  nm,  so 
that  the  correct  size  dependence  of  fracture  could  not  be  predicted. 
On  the  other  hand,  with  high  Gc,  fracture  does  not  occur  even  for 
large  Si  NPs.  With  the  void  radius  of  1.5  nm  and  spacing  of  3  nm, 
cohesive  elements  on  the  surface  completely  lose  their  load¬ 
carrying  capacity  at  a  NP  radius  of  400  nm,  but  not  at  300  nm. 
The  approximated  Gc  is  equal  to  9.4  [J  m-2],  which  shows  good 
agreement  with  recent  experimental  results  [32]. 

4.  Experiments 

The  predicted  size  effect  on  fracture  of  crystalline  Si  NPs  dis¬ 
cussed  above  can  be  compared  with  experiment  using  lithiated  NPs 


with  various  crystalline  axial  orientations  and  sizes.  The  arrays  of  Si 
NPs  were  fabricated  by  dry  etching  of  a  single  crystalline  Si  wafer 
and  subsequently  lithiated;  fracture  was  observed  in  the  scanning 
electron  microscope  (SEM)  as  reported  previously  [15].  For  the 
lithiation,  a  half-cell  of  the  pillar  array  and  a  lithium  foil  was  made 
and  the  voltage  was  swept  down  to  10  mV  (vs.  Li/Li+)  with  a  sweep 
rate  of  0.1  mV  S-1  and  then  held  for  10  h  to  reach  fully  lithiated 
state.  The  tested  axial  crystalline  orientations  of  the  Si  NPs  are 
<100>  and  <110>  and  their  diameters  were  varied  from  140  nm  to 
490  nm.  The  results  for  the  <111  >  pillars  are  taken  from  our  pre¬ 
vious  study  [18].  Fig.  6  shows  representative  SEM  images  of  lithi¬ 
ated  Si  NPs  having  three  axial  orientations  with  both  small  and 
large  diameters  to  reveal  the  size  effect  on  fracture.  The  insets  show 
the  pristine  pillars  with  the  same  magnifications.  All  pillars 
consistently  expand  to  specific  shapes  corresponding  to  their  axial 
orientations  because  of  the  anisotropic  lithiation  behavior  of  crys¬ 
talline  Si,  as  reported  in  our  previous  study  [15,18].  The  relatively 
small  pillars  with  diameters  less  than  290  nm  expand  without 
cracking,  as  shown  in  the  left  hand  column  of  Fig.  6.  In  contrast, 
most  of  the  pillars  with  diameters  larger  than  390  nm  fracture  and 
cracks  propagate  toward  the  center  of  the  pillar,  as  shown  in  the 
right  hand  column  of  Fig.  6. 

For  a  further  investigation  of  the  size  effect  on  the  fracture  of  Si 
NPs  during  lithiation,  the  fraction  of  cracked  pillars  was  counted 
using  SEM.  For  each  of  the  three  axial  orientations  and  for  four 
pillar  sizes,  more  than  100  pillars  were  counted.  Then,  the  “fracture 
ratio”  was  defined  as  dividing  the  number  of  fractured  NPs  by  the 
total  number  of  pillars  counted.  Fig.  7a  shows  the  fracture  ratio  of 
<100>  Si  NPs  of  various  diameters  as  a  column  chart.  When  the 
average  diameters  of  the  pillars  are  200  and  290  nm,  cracks  were 
not  found  in  the  pillars  and  the  overall  fracture  ratio  was  0%  after 
lithiation.  In  contrast,  larger  pillars,  with  390  and  490  nm  di¬ 
ameters,  showed  severe  fracture  (100%  fracture  ratio)  after  lith¬ 
iation.  The  <110>  Si  pillars  exhibit  a  similar  size  effect,  as  shown  in 
Fig.  7b.  When  the  average  diameters  of  the  pillars  are  260  and 
290  nm,  cracking  was  not  found  in  the  pillars.  But  the  pillars  with  a 
diameter  of  360  nm  showed  a  fracture  ratio  of  10%  while  the  pillars 
390  nm  in  diameter  showed  severe  fracture  (100%  fracture  ratio) 
after  lithiation.  The  fracture  behavior  of  <111  >  Si  pillars  shown  in 
our  previous  report  also  exhibits  a  size  effect  18].  For  pillars  with 
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Fig.  6.  SEM  images  showing  the  size  effect  of  fracture  of  Si  nanopillars  after  lithiation  at  10  mV  for  10  h.  Insets  show  corresponding  pristine  pillars  with  same  magnification,  (a) 
290  nm  diameter  <100>  pillar,  (b)  390  nm  diameter  <100>  pillar  <100>  pillars  expand  to  a  cross  shape,  (c)  290  nm  diameter  <110>  pillars,  (d)  390  nm  diameter  <110>  pillar. 
<110>  pillars  expand  to  an  ellipse  shape,  (e)  240  nm  diameter  <111  >  pillar,  (f)  390  nm  diameter  <111  >  pillar.  <111  >  pillars  expand  to  a  circular  shape.  All  scale  bars  are  500  nm. 


diameters  of  140  and  240  nm,  only  a  few  fractures  were  found  after 
lithiation  and  the  overall  fracture  ratios  were  2  and  0%,  respectively. 
The  larger  pillars  with  diameters  of  360  and  390  nm  showed  severe 
fracture  in  most  of  the  pillars  and  their  fracture  ratios  were  88  and 
100%,  respectively.  Based  on  this  study  of  the  fracture  ratios,  it  may 
be  concluded  that  critical  diameter  for  fracture  of  crystalline  Si 
pillars  during  lithiation  is  between  300  and  400  nm,  regardless  of 
the  axial  orientation. 

5.  Discussion 

To  model  the  anisotropic  volume  expansion  in  Si  NPs  during 
lithiation,  various  approaches  have  been  proposed.  One  approach 
involves  using  an  anisotropic  diffusivity  to  mimic  the  anisotropy  of 
the  interface  reaction  rate  [12,19],  coupled  with  the  additional 
assumption  of  an  anisotropic  expansion  coefficient  [33].  Another 
approach,  the  one  we  have  followed  here,  is  based  directly  on  the 
established  anisotropic  phase  boundary  motion  [17].  A  limitation  of 
the  former  approach  is  that  invoking  an  anisotropic  diffusivity  can 
lead  to  a  misunderstanding  of  the  physical  processes  involved, 
since  the  diffusivity  in  a  cubic  structure  is  isotropic.  In  addition, 
amorphization  due  to  lithiation  should  not  have  any  directionality 


so  that  an  anisotropic  expansion  coefficient  is  also  misleading.  For 
these  reasons,  we  have  followed  the  work  of  Pharr  et  al.  17]  and 
directly  modeled  the  anisotropic  expansion  with  anisotropic  phase 
boundary  motion. 

Due  to  the  stress  concentration  that  develops  through  the 
anisotropic  expansion,  the  maximum  stress  would  be  expected  to 
exceed  the  yield  stress  even  without  hardening,  because  of  the 
triaxiality  of  the  stress  state;  the  hydrostatic  stress  does  not  affect 
the  von  Mises  yield  criteria.  However,  the  maximum  stress  should 
be  still  limited  by  the  yield  stress  both  in  anisotropic  and  isotropic 
cases,  because  a  crack  was  observed  to  start  to  form  at  the  surface 
where  the  surface  traction  should  be  free.  Instead,  the  main 
reason  for  easy  fragility  in  anisotropic  expansion  cases  would 
result  from  the  strain  localization  near  the  edge  of  the  phase 
boundaries,  as  shown  in  Fig.  3  and  Supplementary  movies. 
Following  the  suggestion  of  Liang  et  al.  [33  ,  we  conclude  from  this 
that  the  fragility  of  the  lithiated  Si  NPs  stems  from  the  anisotropy 
of  the  lithiation  reaction.  Liang  et  al.  [33]  showed  that  lithiated  Ge 
nanoparticles  are  much  more  fracture  resistant  and  concluded 
that  this  fracture  resistance  derives  from  the  near  isotropy  of  the 
lithiation  reaction.  This  is  also  consistent  with  our  present 
findings. 
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b  Diameter  (nm) 


c  Diameter  (nm) 


Diameter  (nm) 

Fig.  7.  The  effect  of  the  size  on  the  fracture  ratio  for  nanopillars  with  three  kinds  of 
axial  orientations,  <100>,  <110>,  and  <111  >.  (a)  Column  chart  showing  the  fracture 
ratio  for  <100>  nanopillars  of  different  initial  diameters  after  lithiation.  Nanopillars 
with  390  nm  and  490  nm  initial  diameter  have  a  high  fracture  ratio  ( ~  100%),  while 
200  nm  and  290  nm  diameter  pillars  have  zero  fracture  ratio,  (b)  Column  chart 
showing  the  fracture  ratio  for  <110>  nanopillars  of  different  initial  diameters  after 
lithiation.  Nanopillars  with  390  nm  initial  diameter  have  a  high  fracture  ratio  ( ~  100%), 
while  260  nm  and  290  nm  diameter  pillars  have  zero  fracture  ratio.  360  nm  diameter 
pillar  shows  slight  increment  of  fracture  ratio  ( ~  5%).  (c)  Column  chart  showing  the 
fracture  ratio  for  <111  >  nanopillars  of  different  initial  diameters  after  lithiation. 
Nanopillars  with  360  nm  and  390  nm  initial  diameter  have  a  high  fracture  ratio 
(>88%),  while  140  nm  and  240  nm  diameter  pillars  have  low  fracture  ratio  (<~5%). 


For  the  calculation  of  the  critical  size,  in  our  previous  study  34], 
we  were  able  to  compute  the  strain  energy  release  rate,  using  the  J- 
integral  method  [35],  for  pre-existing  sharp  cracks  of  all  possible 
lengths  and  for  all  degrees  of  lithiation.  By  comparing  the 
maximum  possible  strain  energy  release  rate  with  the  work  of 
fracture  of  Si,  we  could  determine  whether  a  NP  of  that  particular 
size  would  fail.  However,  this  kind  of  technique  cannot  be  used  for 
the  present  problem,  because  it  does  not  model  the  crack  propa¬ 
gation.  In  addition,  the  critical  size  based  on  this  technique  is  only  a 
rough  estimate  because  it  is  based  on  a  consideration  of  the  worst 
possible  scenario  for  fracture. 


Instead  of  using  the  J-integral,  cohesive  elements  have  been 
used  to  capture  the  size-dependence  of  fracture,  based  on  the 
microscopic  description  of  void  growth-coalescence.  However, 
recent  in-situ  TEM  experiments  show  that  once  a  crack  is  initiated 
at  the  surface,  it  appears  to  grow  rapidly  to  the  center  [19,20], 
which  behavior  was  not  observed  in  our  finite  cohesive  zone 
modeling.  Instead  of  brittle-like  fracture,  as  seen  in  the  experi¬ 
ments,  the  simulations  show  that  a  crack  does  initiate  at  the  sur¬ 
face,  but  it  does  not  propagate  as  much  as  in  the  experiments.  This 
discrepancy  might  be  understood  by  the  fact  that  newly-formed 
fractured  faces  could  provide  a  new  surface  diffusion  path  for 
lithium  atoms  to  reach  the  pristine  Si,  which  would  alter  the 
morphology  of  the  fracture  process.  Our  fracture  modeling  does  not 
allow  for  such  new  lithium  transport  pathways.  In  addition, 
McDowell  et  al.  [20  reported  that  the  reaction  front  motion  could 
slow  down  due  to  a  stress  effect  on  the  reaction  kinetics  at  the 
phase  boundary.  Newly-formed  fracture  surfaces,  offering  a  new 
path  for  lithium  diffusion,  would  totally  alter  the  stress  state  from 
the  isotropically  well  constrained  stress  state,  so  that  the  observed 
slowing  of  the  reaction  front  would  not  be  observed. 

In  our  microscopic  model,  we  model  a  crack  which  is  initiated 
from  the  surface  and  propagates  to  the  center,  as  observed  in  the 
TEM  observation.  Because  a  crack  occurred  only  in  the  fully  lithi- 
ated  Si,  the  traction-separation  law  would  be  only  for  fully  lithi- 
ated  Li-xSi.  To  explore  the  effect  of  lithium  density  in  the  traction- 
separation  law,  our  model  could  be  extended,  by  making  the 
microscopic  parameters  (void  size  and  spacing)  as  functions  of 
density.  However,  it  often  requires  very  difficult  experimental 
setups  to  get  these  parameters.  In  addition,  recent  experimental 
measurement  showed  that  the  fracture  energy  does  not  vary 
significantly  with  lithium  concentration. 

6.  Conclusions 

We  make  a  phenomenological  model  for  the  anisotropic  lith¬ 
iation  in  crystalline  Si  NPs,  which  is  governed  by  the  kinetics  of 
chemical  reaction  of  lithium  and  pristine  Si.  From  the  model,  similar 
topological  changes  could  be  modeled  and  stress  concentration 
would  occurs  at  the  same  positions  where  fracture  was  observed  in 
experiments.  In  addition,  using  cohesive  elements  associated  with 
microscopic  description  of  void  growth-coalescence,  we  could  es¬ 
timate  work  of  fracture  of  LixSi  alloy,  which  gives  rise  to  the  similar 
critical  diameter  in  experiments  of  <100>  oriented  Si  NPs  during 
lithiation.  These  predictions  are  in  agreement  with  in  experiments, 
which  reported  critical  size  in  the  same  range. 
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Appendix  A.  Elastoplastic  deformation  prior  to  void 
formation 

During  elastic  loading,  the  constitutive  law  can  be  expressed  by 

(7 a  =  2fl£n  +  ^^£xx  +  £yy  +  £zz)  (i  —  x,y,z),  (A.1) 

where  A  is  Lame’s  coefficient  and  /i  is  the  shear  modulus.  Be  sub¬ 
jected  to  uni-axial  deformation,  the  strains  in  the  x—y  plane  are 
zero  : 
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£xx  =  £yy  =  0.  (A.2) 

With  these  relations,  the  stresses  are  related  to  the  axial  strain, 
as  follows  : 


&xx  —  Oxx  —  ^£zz->  &zz  —  (2/x  +  h)£zz  (Elastic  regime).  (A.3) 

This  relation  will  be  used  to  describe  the  first  part  of  the  traction 
separation  law,  when  the  hypothetical  tensile  bars  are  deforming 
elastically. 

For  the  traction  separation  law  during  plastic  deformation  we 
use  the  von-Mises  yield  criterion.  According  to  this  criterion  plastic 
deformation  in  the  bar  starts  at 

0zz  =  (2^  +  (A.4a) 


y 


£zz  — 


2fi  ’ 


(A.4b) 


where  Y  is  the  yield  stress.  Now,  the  total  strains  may  be  divided 
into  elastic  and  plastic  components,  as  follows: 

4°tal  =  4lastic  +  4lastic  (i  =  x,y,z).  (A.5) 

Due  to  the  symmetry  in  the  x—y  plane, 

4-tIC  =  4astlc,  4astlc  =  4*tlc,  and  the  volume  conservation 

during  the  plastic  deformation  £^astlc  +  £y]?stlc  +  £zzastlc  =  0,  plastic 
strains  in  the  x—y  plane  can  be  expressed  by 

plastic  _  plastic  _  1  plastic  ^ 

£xx  —  £yy  —  2£zz  ' 

With  the  condition  that  total  strains  in  x—y  plane  are  zero, 
elastic  strains  in  x—y  plane  are  then  expressed  by 


elastic  _  elastic  _  1  plastic 
bxx  ~  byy  ~  2  zz 

Following  Hooke’s  law,  the  stresses  are  expressed  by 


(A.7) 


an  =  2nef^  +  A(4astic  +  4astic  +  4astic)  (i  =  x,y,z). 

(A.8) 

Since  the  yield  condition  is  always  satisfied  during  plastic 
deformation,  the  von-Mises  yield  criterion  leads  to 


plastic  _  9  elastic 

bzz  ^bzz 


Y 

M' 


(A.9) 


From  this  relation,  the  stress-strain  relation  in  the  loading  di¬ 
rection  after  yielding  but  prior  to  void  formation  is  expressed  by 


ozz  =  0^  +  4ztal  + 1  y  (Plastic  regime) . 


(A.10) 


Appendix  B.  Damage  description  with  cohesive  elements 

In  this  work,  we  simulate  the  crack  initiation  and  propagation 
using  cohesive  elements.  Initial  response  of  cohesive  elements  is 
assumed  to  be  linear-elastic  prior  to  satisfying  the  damage  initia¬ 
tion  criterion,  which  is  traction  approached  to  the  specified  stress 
(gd).  Beyond  this  point,  cohesive  elements  start  to  lose  load¬ 
carrying  capacity  by  degrading  their  element  stiffness  and 
completely  fail  and  should  be  deleted  from  the  model  when  the 
separation  distance  reaches  to  Wf.  The  area  beneath  the  traction- 


separation  law  is  equal  to  the  work  of  fracture  (Gc)  required  for  a 
crack  to  propagate  in  the  context  of  the  linear  elastic  fracture  me¬ 
chanics.  It  should  be  noted  that  cohesive  elements  only  respond 
linear-elastically  under  compression,  because  fracture  only  occur 
under  tension  or  shear  loading. 
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List  of  symbols 


a  [nm]:  radius  of  the  Si  core  in  a  pressurized  tube 
b  [nm]:  outer  radius  of  a  pressurized  tube 
r  [nm]:  radial  position  in  a  pressurized  tube 
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Y  [nm]:  yield  stress 

u  [nm]:  radial  displacement 

w  [nm]:  separation  displacement  in  cohesive  zone 

wP  [nm]:  separation  displacement  at  onset  of  the  plasticity 

wD  [nm]:  separation  displacement  at  void  formation 

Wf[nm]:  maximum  separation  displacement 

Gc  [J  m~2]:  work  of  fracture 

Greek 


a  [nm]:  void  radius 
ao  [nm]:  initial  void  radius 
a  [nm]:  rate  of  change  of  the  void  radius 
/?  [nm]:  void  spacing 


X  [GPa]:  Lame’s  coefficient 
/i  [GPa]:  shear  modulus 
y s  [J  m  2]:  surface  energy 
a  [GPa]:  traction  in  cohesive  zone 
up  [GPa]:  traction  at  onset  of  the  plasticity 
(td  [GPa]:  traction  at  void  formation 
dm  [GPa]:  mean  stress 

o'**  ayy,  azz  [GPa]:  stresses  in  Cartesian  coordinate  system 

am  add,  a w  [GPa]:  stresses  in  cylindrical/spherical  coordinate  system 

£eiastic '  Mastic,  £totai  j .y.  strains  (elastic/plastic/total) 

ee,  ep,  el  [s~J]:  strain  rates  (elastic/plastic/transformation) 
i  [s~:]:  effective  strain  rate 

err ,  £dd->  Is  1 1:  strain  rates  in  spherical  coordinate  system 


